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ABSTRACT 


This  report  discusses  various  techniques  of  designing 
auto-adaptive  (self-changing)  filters,  and  the  application 
of  one  particular  technique  to  two  teleseisms.  This  report 
also  discusses  the  calculation  of  running  correlation  func¬ 
tions,  using  a  moving  time  window,  and  iterative  solutions 
to  the  multichannel  normal  equations.  A  listing  is  included 
of  a  program  which  will  calculate  an  auto-adaptive  multi¬ 
channel  filter  by  the  use  of  running  correlation  functions 
and  the  method  of  steepest  descent. 


: 


SECTION  I.  INTRODUCTION 


In  the  analysis  of  seismic  records,  it  is  apparent  to  the 
analyst  that  there  are  both  long-term  and  short-term  variations 
in  the  noise  statistics  of  seismic  array  data.  An  example  of 
long-term  variation  in  noise  statistics  is  the  seasonal  change 
in  the  noise  power  level;  an  example  of  short-term  variation 
is  the  short  bursts  of  obviously  coherent  energy  which  appear 
on  most  seismic  records  which  the  analyst  sees.  The  question 
arises  as  to  whether  significant  rejection  of  time-varying 
noise  can  be  achieved  by  multichannel  filtering. 

A  filter  whose  coefficients  depend  upon  the  time-varying 
statistics  of  the  nois<n  is  termed  an  adaptive  filter.  A  filter 
whicn,  as  the  noise  statistics  change,  automatically  updates 
itself  without  the  need  of  human  intervention  or  supervision  is 
termed  an  auto-adaptive  (self- changing)  filter.  The  auto- 
adaptive  filters  described  in  this  report  were  designed  to  be 
linear  but  time-varying  multichannel  filters. 


SECTION  II.  THEORY  OF  MULTICHANNEL 
ADAPTIVE  WIENER  FILTERING 

1-  The  Mean  Square  Error  Criterion 

The  multichannel  wave  shaping  problem  as  illustrated  in 
Figure  1  may  be  stated  as  follows: 


Figure  1.  The  multichannel  wave  shaping  problem 


Given  a  set  of  n  input  time-series  data  traces  x^  through  x  , 
and  a  set  of  m  desired  output  traces  d1  through  dm,  design  a  set  of 
n  filters  to  transform  the  n  input  data  series  into  a  set  of  approx 
imations  y.±  through  ym  of  the  m  output  data  series,  subject  to  the 
constraint  that  the  mean  square  error  (the  time  average  of 

?  2 

Z  (d,*  “  y,* )  )  be  a  minimum. 
i=l  1  1 


folloT.  Pr0bl8m  “ay  be  SUCCin0tly  State*  in  notation  as 

The  multichannel  convolution  matrix  (X)  of  the  set  of  input 

s:ipled  by  the  fut~ — <»  ^  p 

In  vector  notation 
X  F  =  Y 

2.1 

where  X  is  the  multichannel  analog  of  the  sincri*  i 

truiar*  g  single  channel  rectan¬ 

gular  Toeplitz  convolution  matrix. 

error  vector  equals  the  difference  between  the  desired 
output  vector  and  the  actual  output  vector: 

£  =  D  -  Y 

2.2 

The  mean  square  error  is  the  expected  value  of  £l£.  The  t  super- 
8  -rj.pt  denotes  matrix  transposition. 


E  =  MSE  =  E  (D1  ~  fV)  (D  -  XF)} 

-W 

(2.3)  can  be  expanded  to  give 

E  {£2}  =  EiD^}  +  E  {-fVd  -  dV)  +  E  {fVxF} 

Since  the  transpose  of  a  scalar  is  equal  to  itself, 

(FtXtD)  =  D*XF 
and  (2.4)  can  be  written 

E  ^-E(D}-2E  {F*X*D}  +  E  (F^X^XF)  2 

In  terms  of  R  =  (E  {X**}),  the  multichannel  correlation 
matrix  of  the  input  data,  and  G  i  (E  {X^D} ) ,  the  multichannel 
crosscorrelation  vector  between  the  input  with  the  desired 
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output,  we  have 

E  U2}  =  DJ  -  2  FtG  +  FtRF  5-7 

The  criterion  that  the  mean  square  error  be  a  minimum  leads 
to  the  equation 


=  -  2G  ♦  2RF  2* 8 

TP 

Setting  =  0  gives  the  multichannel  normal  equations : 

aF* 

RF  =  G 

where  the  expected  value  of  the  gradient  of  the  mean  square  error  is 


-  =  -  2  (RF-G ) 

For  the  case  N  =  2,M=  2,  L=  3,  the  multichannel  normal 
equations  take  the  form 


2.10 


2.11 


fn<0>  h?*0’ 


fjj(°)  r?7<o> 


8u<0) 

*7!<01  *??<0> 


fnd)  f17(i) 


f7l(l)  f„(7) 


«uo>  «17C1) 


•7iU>  g77ci) 


fu<?>  f17<?> 


r..(7)  f..(7) 


«n(7)  *,,(7) 


■..(7)  *„ .(7) 


where 


N  =  number  of  input  channels 
M  =  number  of  output  channels 

L  =  n,'mber  of  filter  points  (number  of  lags  in  the  con  elation 
functions ) 

It  should  be  noted  that  the  multichannel  correlation  matrix, 
R,  in  Equation  (2.11)  is  a  symmetric  matrix,  and  moreover  that 
writing  it  in  terms  of  submatrices , 


we  see  that  it  is  a  block  Toeplitz  matrix  whose  cross-diagonal 
terms  are  transposes  (i.e.,  r?.  =  r .  . ) .  This  characteristic 

of  the  R  matrix  leads  to  some  simplification  in  the  solution 
of  the  equations . 


The  principle  of  orthogonality  (Papoulis,  1965)  states 
that  lor  the  optimum  multichannel  filter,  the  expected  value  of 
the  convolution  of  the  error  trace  and  the  input  data  channels 
is  zero  (for  X  or  E  having  zero  mean),  i.e., 


E  {X1^}  =  0 


2.13 


We  have  shown  (Equation( 2 . 9 ) )  that  the  optimum  wave 
shaping  filter  satisfies  the  equation 

(RF  -  G)  =  0 
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4 


Equation  (2.14)  can  be  factored: 

E  (Xt  XF  -  XtD)  =  0  2.15 

or 

E  (Xt  (D  -  XF) }  =  0  2.16 

or 

E  {Xt(5)>  =  0  2.17 

hence 

E  (X*  £}  =  G-RF  2.18 

*0 

Equation  (2.18)  and  (2.14)  ere  thus  equivalent.  Therefore, 
the  gradient  of  the  mean  square  error,  -2(RF-G),  is  also  given 
by  2E  {X1^}.  Widrow  (1966)  makes  use  of  this  equivalence  to  de¬ 
sign  adaptive  filters.  As  a  function  of  the  iteration  number  n, 
the  "method  of  steepest  descent,"  or  point-slope  iterative  pro¬ 
cedure  for  the  solution  of  the  Wiener-Hopf  equation  in  the  time 
domain,  may  be  written 

Fn+1  =  Fn  -  k  E  (xH)  2.19 

or 

pnvl  _  pn  _k(RFn  _  q)  2.20 

where  k  is  the  magnitude  of  the  step  taken  in  the  direction 
opposite  to  the  gradient  of  the  rean  square  error. 

Widrow  does,  however,  make  ont.  further  assumption,  namely 
that  the  expected  value  of  the  convolution  of  the  input  data 
channels  with  the  error  trace  is  equal  to  the  product  of  the 
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input  with  the  error  trace  at  any  one  instant  of  time.  This 
assumption  leads  to  a  very  noisy,  but  unbiased,  estimate  of 
the  gradient  of  the  mean  square  error. 


3.  Adaptive  Filtering  Solutions 

The  solution  to  the  multichannel  Wiener-Hopf  equation 
RF  =  6  may  be  obtained  by  one  of  three  methods : 


1)  direct  inversion  of  the  R  matrix; 

2)  application  of  the  Levinson-Robinson  recursion  formula 
(Levinson,  1949;  Robinson,  1967); 

3)  application  of  an  iterative  technique. 


The  method  of  steepest  descent  (Equations (2, 19)  and  (2.20)  as 
well  as  other  iterative  solutions,  depends  on  the  fact  that  the 
mean-square  error  function  is  a  quadratic  function  (Equation  2.6) 
of  the  filter  weights.  The  mean-square  error  function  has  a 
unique  minimum.  Figure  2  illustrates  the  quadratic  mean- square 
error  surface  as  a  function  of  a  two  point  filter.  The  pro¬ 
jection  of  lines  of  constant  mean-square  error  on  the  MSE  surface 
are  conic  sections  in  the  filter  space.  The  gradient  of  the 
scalar  function  in  filter  space  is  -  .  As  the  noise  statistics 

change,  one  can  picture  the  mean-squa&e-error  surface  translating, 
rotating,  changing  its  ellipticity,  opening  and  closing. 


The  taking  of  a  step  in  the  opposite  direction  of  the  grad¬ 
ient  in  filter  space  is  equivalent  to  using  Equation  (2.19)  or 
(2.20). 


There  are  three  other  iterative  solutions  (Varga,  1965)  under 
investigation  at  the  Seismic  Data  Laboratory, 

a)  Iterative  solution  No.  1  is  derived  as  follows: 


RF  =  G 


2.21 


(R+CI-CI )F  =  G 

2.22 

CIF  =  G-CR-CI )F 

2.23 

F  =  C”1(G-RF )  +  F 

2.24 

70+1  =  C"1(G-RFm)  ♦  Fm 

2.25 

■m+l  _  pm  +  c“1(G-RFm) 

2.  26 

where  C  is  arbitrary.  However,  C  should  probably  be  rQ  in 

order  that  the  filter  should  not  be  changed  by  values  greater 
than  unity. 

It  should  be  pointed  out  that  this  iterative  solution  is 
equivalent  to  Equation  (2.20)  since  RF-G  is  the  gradient  of  the 
mean  square  error . 

b)  Iterative  solution  No.  2,  is  a  point- Jacobi  technique  in 
which  R  is  split  into  a  diagonal  matrix,  an  upper  triangular  matrix, 
and  a  lower  triangular  matrix: 


R  =  (rQI  +  U  +  L) 


RF  =  G 


(rQI  +  J  +  L)F  =  G 


rQF  =  G  -  (U+L)F 
Fm+1  .  r  "1{G  .  (U+L)F“) 


This  method  is  quite  simple  to  implement  since  can  be 

easily  stored  in  the  computer. 


2.  27 

2.28 

2.29 

2.30 

2.31 


•  -i 


( 


' 

•  ’ 


1  I 


'  *  M 

•  'm : 


*V  * 

9  f  » 


ft 


c)  Iterative  solution  No.  3  is  the  Gauss-Sf;idel  method: 


R  =  (rQI  +  U  +  L) 


2.32 


RF  =  G 


2.33 


<rQI  +  U)F  =  G  -  LF 


2.34 


,n+1  =  (rQI  +  U)”1  (G  -  LFn) 


2.35 


The  inverse  of  an  upper  or  lower  Toeplitz  triangular  matrix  was 
shown  by  Meserve  (1968)  to  be  an  easily  obtainable  upper  or  lower 
Toeplitz  triangular  matrix.  The  derivation  for  a  4  x  4  case  is: 


t  t  t 

r0  rl  r2  r3 


xQ  xx  *2 


10  0  0 


0  r0  rl  r2  x  0  xo  X1  x2 


0  I  0  0  2.  36 


r0  rl 


x0  X1 


0  0  10 


0  0  r, 


0  x, 


0  0  0  1 


then 


roxo  =  11  xo  =  ro 


2.37 


rl  xo  +  roxi  =  °>  X1  1  -V1  ri  xo 


'x0  rl  x0 


2.38 


t  t  t  t  t 

TqX2  +  r£x^  +  r2x0  =  x^  =  +  xQrixorixO”xOr2xO  etc. 


2.  39 


Thus  we  have  a  simple  recursive  scheme  for  finding  the  inverse 
of  any  upper  (or  lower)  Toeplitz  triangular  matrix. 

We  constructed  a  test  case  for  the  Gauss-Seidel  method  where 
the  number  of  input  channels  was  2,  the  number  of  lags  was  3  and 
the  number  of  outputs  was  2.  The  equation  to  be  solved  is 


t  t 

" 

r0  rl  r2 

f0 

g0 

rl  ro 

X 

fl 

s 

gl 

_r2  rl  r0 

f  2 

_g2 

where 


2.41 


S 
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The  initial  guess  for  the  filters  was 


The  results  for  this  test  case  are  shown  in  Figure  3.  The 

plotted  points  are  the  filter  points  corresponding  the  second 
column  of  Equation  (2.42). 

The  Gauss-Seidel  method  is  guaranteed  to  converge  for 
diagonally  dominant  matrices  (Varga,  1965).  By  diagonally 
dominant,  we  mean  that 


r. . (k) 
ID 


<Cr..(0)  r..(0)]1/2; 


for 


The  stipulation  that  a  matrix  be  diagonally  dominant  is  no 
real  restriction  in  the  analysis  of  time-series  data  since  diagonal 
dominance  implies  that  no  input  channel  can  be  more  highly  corre- 
lated  with  another  channel  than  with  itself. 

In  order  to  use  any  of  the  three  iterative  solutions  (which 
depend  upon  knowing  R  and  G)  in  an  adaptive  filter  program  it  is 
necessary  to  have  some  means  of  updating  the  multichannel  correla¬ 
tion  matncies  with  time.  A  method  of  updating  the  R  matrix 
and  the  G  matrix  with  time  is  presented  below. 


It  is  possible  to  update  the  correlation  matrix  of  a  multi¬ 
channel  time  series  with  the  arrival  of  each  new  data  point. 

Given  a  multichannel  time  series  X.(t).  i  =  l  9  « 

‘  )  •  .  .  , n. 
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O&Hhsds  XJDWtn 


ro 

o 

is 

•H 

Cn 


^  H  ,4  H  W  OS  UOWfcHUHWJJHW 


t=t ' 


t=t ' +LF 


Consider  the  possibility  of  taking  all  possible  correlations  of 
all  the  data  points  in  the  time  window  of  length  LF  (LF  =  "leng 
of  fitting  interval"). 


r 


t'+LF 
•  • 


(k) 


t'+LF 

l  X . (s )X. (s+k) 
s=t'  1  3 


2 


Consider  what  happens  to  the  correlation  matrix  when  the  time 
window  is  advanced  one  point. 


tLF+1(k> 


t'+LF+l 

l  X. (s)X. (s+k) 

s=t'+l  1  3 


2, 


rt'+LF+l(k)  .  rt'+LF(k)  _  x.(t' )X.(t'+k)  +  X. (t'+LF+l-k)* 
i]  13  ID  i 

Xj  (t=+LF+l) 


2.48 


Equation  (2.48)  states  that  the  ij'th  element  of  r(k)  after 
a  new  data  p'jint  arrives  is  equal  to  the  previously  calculated 
element  minu-3  one  old  term  plus  one  new  term. 

Note  that  in  the  last  term  of  the  preceding  expression,  the 
lag  parameter  k  was  shifted  to  X^  from  X^  and  the  sign  of  k  was 
changed.  This  was  done  since  data  points  ahead  of  t  =  t'+LF+l 
are  as  yet  unknown. 

In  the  solution  of  the  Wiener- Hopf  equation  RF  =  G,  R  is 
the  known  correlation  matrix,  F  is  the  unknown  filter  vector 
and  G  is  the  crosscorrelation  vector  of  the  desired  output  with 
the  input  channels.  We  can  update  an  element  of  G  as  easily  as 
an  element  of  R  by  simply  replacing  X^  in  the  Equations  (1.1)  - 
(1.4)  by  d^  and  replacing  r  by  g: 


Hence 


gt^+LF+l(k)  .  gt^+LF(k)  _  di(t»)  x. (t’+k)  +  d^t'+LF+l-k)- 


2.49 


Xj  (t'+LF+l) 


In  summary,  there  appear  to  be  two  basic  methods  of  linear 
multichannel  time-varying  filtering. 

Method  1  is  to: 


a)  solve  the  Wiener-Hopf  equation  once  to  get  the  exact 

solution  for  F  (or  alternatively  start  with  an  arbitrary 
guess ) ; 
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b)  update  R  and  G  by  the  method  in  Equations  (2.48)  and 
(2.49); 

c)  periodically  perform  one  of  the  iterative  techniques  to 
update  the  filter;  or  else  use  some  parameter  such  as  the 
variation  of  the  mean  square  error  to  determine  the 
necessity  to  update  the  filter. 

Method  2  depends  upon  the  estimation  of  the  gradient  of  the 
mean  square  error: 

a)  start  with  an  arbitrary  set  of  filter  coefficients 

b)  estimate  RF  -  G,  or  equivalently,  estimate  E{XtC) 

c)  then  Fra+1  =  Fm  -  k  (RFm-G )  or  Fm+1  =  Fm  -  k  EfX1^) 

The  technique  used  by  Widrow  (1966)  and  Burg  et:  al  (1967) 
is  to  estimate  the  expected  value  E{Xt^}  by  X£.  The  technique 
used  in  this  report  estimates  E{X£)  by  use  of  a  running  cross- 
correlation  vector  between  the  input  and  the  error  trace,  updated 
with  the  arrival  of  each  new  data  point  by  the  method  of  Equation 
(2.49). 
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SECTION  III.  EXPERIMENTAL  RESULTS 


In  order  to  use  the  steepest  descent  algorithm  Equation  (2.19), 
it  is  first  necessary  to  get  some  desired  output  trace  from  the 
set  of  input  data  channels.  An  extremely  useful  multichannel  fil¬ 
ter  is  the  maximum-likelihood  filter  (Kelly  and  Levin,  1964)  which 
attempts  to  pass  signals  which  are  the  same  on  every  channel  and 
to  reject  signals  (or  noise)  which  are  not  the  same  on  every  channel 
(For  a  derivation  of  the  maximum  likelihood  filter  see  McCowan, 
1968).  The  maximum-likelihood  filter  may  be  cast  into  terms  of 
a  prediction-error  filter  (Burg  et  ad,  1967),  For  the  sake  of 
clarity  and  coherence,  the  algorithm  designed  by  Burg  et  al  (1967) 
is  presented  here. 

The  maximum  likelihood  filter  is  a  multichannel  filter  that 
is  designed  subject  to  the  constraint: 


N 

1  f-« 

i  =  l  13 


6. 

3  8 


3.1 


or 


N 


f,.  =  5 


.  -  l 

3 s  it 


i=  2 


V 


3.2 


Multichannel  convolution  of  the  filter  with  the  inputs  gives 
the  output: 


Vs 


N  L 

i=l  k=l  1JC  13  * 


or 


Vs 


L 

l  f 


lk  Xlj-k 


L  N 

+  y  y  f ..  x. .  . 

L  -L  -  ik  lj-k 


k=l  i=  2 


Substituting  (3.2)  into  (3.4)  gives 


Vs 


L  f  N  L  N 

Jl  ps  "  il2  fikJ  Xlj-k  +  Jl  i=  2  VV*’ 


3.3 


3.4 


3.5 
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3.6 


Y. 

]-s 


Equation  3.6  is  seen  to  be  the  error  remaining  in  the  attempt 
to  predict  the  (j-s)’th  value  of  from  che  differences  between 
X1  and  each  of  the  rest  of  the  channels.  The  filters  designed  in 
this  section  have  made  use  of  this  algorithm.  The  cross¬ 
correlation  between  the  input  and  the  error  trace  was  calculated 
and  updated  with  the  arrival  of  each  new  data  point  by  the  method 
of  Equation  (2.49).  The  new  filter  for  each  new  data  point  was 
calculated  from  Equation  (2.19).  The  original  guess  for  the 
optimum  filter  was,  in  both  cases,  1/N,  0  ,  0,  0...  for  each  of 
the  N  channels  used  to  predict  the  (N  **  l)’st  channel.  The  data 
was  detrended,  demagnified,  time-aligned  and  bandpassed  filtered 
by  the  "SDL  Bandpass  Filter"  (Figure  4).  The  length  of  the 
adaptive  filter  was  30  points  (1.5  sec).  The  best  results  were 
obtained  by  shifting  the  trace  to  be  predicted  20  points  (1  second 
ahead  in  time).  This  time  shift  corresponds  to  designing  a  fil¬ 
ter  to  predict  20  points  ahead.  The  filters  were  one-sided 
filters  which  operate  only  on  past  and  present  data.  The  length 
of  the  fitting  interval  was  320  points  (16  seconds).  The  auto- 
adaptive  filters  in  this  report  were  capable  of  operating  in  one 
of  three  modes,  namely: 

Mode  1:  Operate  (filter  the  data),  learn  (update  the 
filter  with  each  new  data  point),  and  update  (update  the 
gradient  of  the  mean  square  error  with  each  new  data 
point ) . 


Mode  2:  Operate  and  update  only. 


Mode  3:  Operate  only. 

The  data  were  stored  in  core  in  a  block-multiplexed  form 
(512  points  from  each  channel).  Thus  there  were  192  points 
outside  the  fitting  interval  which  could  be  operated  on  without 
reading  new  data  points  from  the  disc.  The  filters  were  designed 
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to  operate  an  Mode  1  for  100  points  and  in  Mode  2  for  92  points, 
alternatively.  After  the  signal  arrived  the  filters  were  con¬ 
strained  to  operate  in  Mode  2  for  all  192  points.  In  an  online 
si  uation  it  would  be  necessary  to  have  in  conjunction  with  the 
a  aptive  filter  a  signal  detector  system  which,  after  a  signal 
was  detected,  would  cause  the  adaptive  filter  to  operate  only  in 
ode  until  a  preset  time  after  the  signal  arrival.  However, 
even  if  a  signal  is  not  detected,  the  gradient  of  the  mean  square 
error  win  not  change  drastically  until  the  signal  is  well  within 
the  fitting  interval,  since  the  scale  factor  k  in  Equation  (2.20) 
was  made  to  depend  upon  the  zero  lag  of  the  autocorrelation  of 
each  trace.  The  zero  lag  of  the  autocorrelation  of  each  trace 

was  updated  with  the  arrival  of  each  new  data  point  in  a  manner 
similar  to  that  of  the  gradient: 


3.  7 


3.8 


rT.+LF+1(0)  =  r^  +  LF(0)  -  X*  Ct')  +  X?  (t'+LF+l) 

and  thus  the  scale  factor  k  in  Equation  (2.19)  for  the  i-th 
channel  was : 

k(i)  =  k/r*!+LF+1(0) 

11 

In  other  words,  the  magnitude  of  the  step  taken  opposite  the 
irection  of  the  gradient  of  the  mean  square  error  was  inversly 
proportional  to  the  magnitude  of  the  zero  lag  of  the  auto¬ 
correlation  of  each  channel.  This  procedure  was  equivalent  to 

normalizing  the  variance  of  each  channel  to  unity  over  the  moving 
time  window.  & 

Figure  5  is  a  map  of  the  ten  sensors  taken  from  the  LASA 
subarray.  Figure  $  is  the  array  response  of  the  partial  D1 
subarray.  The  data  and  the  array  used  are  the  same  used  in  a 
previous  report  (McCowan,  1988),  in  order  to  compare  directly  the 
results  of  adaptive  filtering  with  time  stationary  filtering. 
Table  III-2  is  a  list  of  the  sensor  positions.  Table  III- 3  is  the 
pertainent  data  for  the  Aleutian  event.  Table  III- 4  gives  the 
contour  increment  for  the  f-k  response  function. 
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Figure  6.  Array  impulse  response  function 
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TABLE  III-2 


MULTICHANNEL  FILTER  REPORT 

Array  configuration 

Channel 

No. 

Channel 

ID 

X  Coordinate 
(km) 

Y  Coordinate 
(km) 

Demagnification 
(counts/ 10” 4 
microns) 

1 

10 

0.000 

0.000 

2.55 

2 

81 

-0.716 

3.426 

3.09 

3 

32 

0.745 

0.667 

2.78 

4 

52 

1.610 

1.440 

3.18 

5 

83 

3.325 

-1.093 

3.07 

6 

34 

0.204 

-0.979 

2.48 

7 

54 

0.409 

-1.958 

3.00 

8 

85 

-2.609 

-2.333 

3.05 

9 

36 

-0.950 

0.312 

2.69 

10 

56 

-1.900 

0.625 

2.75 

TABLE  II I- 3 


Location  of  Event 


Date : 
Region: 
Magnitude* : 


10  November  1965 
Aleutian  Islands 
4.1 


Origin  Time*:  03:48:78.0  GMT 
Latitude*:  51.2°  N 

Longitude*:  178.8°  W 

Focal  Depth*:  "33  km" 

Epicentral  Distance  to  LASA  A010: 


Azimuth  of  Epicenter  From  LASA  A010: 
Surface  Velocity  at  LASA  A010: 


46.1°  (5126  km) 
304° 

14  km/ sec 


*Data  from  USCSGS  Preliminary  Determination  of 


Epicenter. 


Figure  7  is  the  result  of  applying  the  auto-adaptive  filter 
to  the  input  data.  The  first  trace  is  a  timing  trace  with 
1- second  pips,  every  fifth  one  of  which  is  accentuated.  Traces 
2  through  11  are  the  observed  data,  trace  11  being  the  one  to 
be  predicted;  this  trace  is  shifted  1  second  to  the  left.  Trace 
12  is  the  output  of  the  auto-adaptive  filter.  Trace  13  (the 
error  trace)  is  the  difference  between  trace  11  and  trace  12. 
Trace  14  is  the  sum  of  the  nine  input  traces  used  to  predict 
trace  11.  The  output  of  the  adaptive  filter  is  seen  to  be 
identical  to  the  sum  at  the  start  of  the  record  (as  it  should 
be,  since  the  initial  guess  for  the  filters  was  1/N,  0,  0,  0...). 

Figure  8  is  a  continuation  in  time  of  the  same  data.  We 
can  see  that  the  filter  is  indeed  learning  from  the  input  data, 
since  the  noise  power  is  getting  progressively  less  in  trace  12 
as  compared  to  trace  14.  In  general  the  same  peaks  seem  to  be 
coming  through,  but  the  power  is  quite  a  bit  less  in  the  output 
of  the  adaptive  filter  (trace  12)  as  compared  to  the  beamed  sum 
trace  (trace  14)  as  time  goes  on.  Figure  9  is  a  continuation 
of  the  same  data.  The  first  motion  of  the  signal  seems  to  be 
enhanced  as  does  a  secondly  phase,  probably  pp.  The  signal- 
to-noise  ratio  improvement  l  2  db  over  the  bandpass  phased 
sum,  due  to  some  signal  rejection  by  the  adaptive  filter.  The 
second  phase  occurs  roughly  nine  seconds  after  the  primary  phase, 
the  time  delay  being  about  right  for  pP. 

It  should  be  noted  that  the  progressive  decrease  in  the 
noise  power  with  time  is  due  to  two  factors:  one  is  the  con¬ 
vergence  of  the  adaptive  filter  from  the  initial  guess  toward 
the  optimum  time-stationary  filter;  the  other  factor  is  the 
improvement  achieved  by  reducing  time- varying  coherent  noise. 

The  small  improvement  (2  db)  above  the  band  limited  phaued  sum 
seems  to  indicate  that  the  amount  of  time-varying  coherent 
noise  power  is  low.  Of  course,  another  factor  to  consider  is 
the  poor  array  response  (Figure  6)  which  causes  coherent  noise 
to  alias  back  into  the  main  lobe. 

Table  III-5  gives  the  pertainent  data  for  the  W. Andreanof 
Islands  event. 
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Figure  9.  The  Aleutian  Event 


TABLE  1 1 1  -  5 


Location  of  Event 


Date : 
Region*: 
Magnitude*: 
Origin  Time* 
Latitude* : 
Longitude*: 
Focal  Depth* 


22  February  1968 

West  Andreanof  Islands 

3.8 

14:43: 46 
51.  6N 
176. 2W 


6t>  km 

Epicentral  Distance  to  UBQ:  46.137°  (5130  km) 
Azimuth  of  Epicenter  to  UBO:  76.773° 

Surface  Velocity  at  UBO:  14  km/sec 


*  Data  from  USCSGS  Preliminary  Determination  of  Epicenter, 


Figure  10  is  a  map  of  UBO  nine-channel  surface  array. 

Figure  11  is  the  array  response  function.  Figure  12  is  the 
input  data  after  it  has  been  demagnified,  bandpass  filtered, 
and  beamed;  the  various  outputs  are  also  shown.  Trace  1  is  the 
timing  trace  (one  second  pips,  with  every  fifth  pip  accentuated). 
Traces  2  through  11  are  the  input  data.  Trace  11  is  the  pre¬ 
dicted  trace  (shifted  1  second  ahead  in  time).  Trace  12  is  the 
output  of  the  adaptive  filter.  Trace  13  is  the  beam-formed  sum 
and  trace  14  (the  error  trace)  is  the  difference  between  trace 
12  and  trace  13. 

Once  again,  we  note  that  at  the  beginning  of  the  record, 
trace  12  and  trace  13  are  identical  due  to  our  original  estimate 
of  the  optimum  filter  (the  original  estimate  being  1/N, 0,0,0... 
for  each  channel). 

Figure  13  shows  a  continuation  in  time  of  the  same  data. 

Note  that  the  auto-adaptive  filter  is  performing  better  than  the 
phased  sum  as  the  filter  learns  more  about  the  noise  character¬ 
istics.  Figure  14  shows  a  continuation  of  the  same  data.  The 
filter  was  constrained  to  operate  in  Mode  2  shortly  before  the 
expected  arrival  time  of  the  signal  (14; 52: 04).  Although  it 
was  not  practical  to  measure  signal-to-noise  ratios  with  such 
a  small  signal,  it  is  obvious  to  the  eye  that  the  auto-adaptive 
filter  is  doing  a  very  good  job.  In  order  to  compare  the 
results  of  the  auto-adaptive  filter  with  the  on-line  analog 
MCF’s  we  have  reproduced  a  section  of  film  from  Develocorder 
Number  5  at  UBO  (Figure  15).  A  description  of  the  design  of 
various  filters  may  be  found  in  a  report  by  Edwards  (1965). 

Table  III-l  is  a  description  from  a  report  (Leichliter,  1967) 
aboqt  the  on-line- f ilters .  The  auto-adaptive  filter  (Figure  14) 
seems  (by  eye)  to  be  doing  a  better  job  than  MCF  11  (trace  2) 
which  is  an  infinite  velocity  measured-noise  filter.  Of  course, 
it  is  net  quite  fair  to  compare  an  analog  filter  operating  on 
analog  data  with  a  digital  filter  operating  on  digital  data,  since 
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Figure  11.  Array  impulse  response  function 
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all  things  being  equal,  the  digital  filter  will  always  out-perform 
the  analog  filter.  The  interesting  thing  about  Figure  14  and  15 
is  that  the  analog  filters  operating  on  the  vertical  array  seem 
to  be  cut-performing  (in  terms  of  first  motion  enhancement)  even 
the  digital  auto-adaptive  filter.  The  reason  for  the  better  per¬ 
formance  of  the  vertical  array  stems  from  the  fact  that  the  noise 
level  is  lower  on  the  vertical  array,  and  the  fact  that  the 
surface  array  is  small. 

Figure  16  illustrates  the  application  of  a  time— stationary 
maximum  likelihood  filter  to  the  input  data.  Trace  1  is  the  sum 
of  the  ten  input  traces  (beamed  and  bandpassed).  Trace  2  is 
the  output  of  the  maximum  likelihood  filter.  This  maximum  likeli¬ 
hood  filter  was  designed  on  3800  points  of  noise  before  the 
signal.  The  filter  was  a  one  sided  21  point  filter  which  operated 
on  every  other  data  point,  hence  the  filter  was  effectively  a 
2.1  second  filter. 
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SECTION  IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


From  this  report  it  can  be  concluded  that: 

.  The  method  of  updating  the  correlation  functions  used 
ir.  this  report  is  a  valid  method  for  calculating  relatively 
stable  and  slowly  time  varying  correlation  functions. 

.  The  auto-adaptive  filters  designed  from  the  time-varying 
noise  correlation  functions  do  converge,  and  they  yield 
better  results  (2  db  better  on  the  LASA  event)  than  simple 
bandpass  filtering  and  beamforming.  The  concept  of  auto- 
adaptive  filtering  could  be  applied  to  single  channel 
(prediction,  prediction-error,  and  deghost)  filtering  as 
well  as  to  multichannel  filtering. 

From  this  study,  it  can  be  recommended  that: 

.  Additional  iterative  techniques  such  as  the  Gauss-Seidel 
method  be  investigated. 

.  A  single-channel  deghosting  auto-adaptive  filter  program 
be  written  for  vertical  array  data. 

.  A  single  channel  prediction  -  error  auto-adaptive  filter 
program  be  written  to  process  the  output  of  isotropic  multi 
channel  time  stationary  filters  or  to  process  beamed  sum 
traces . 
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APPENDIX:  PROGRAM  LISTINGS 

This  section  of  the  report  presents  the  program  listing  for 
Program  Adamax3  and  Subroutine  Upgrad .  Program  Adamax3  performs 
the  operations  for  adaptive  filtering  in  Mode  1  and  Mode  2.  Sub¬ 
routine  Upgrad  performs  the  actual  updating  of  the  gradient  of 
the  mean  square  error  and  the  actual  filtering.  Both  routines 
are  written  for  clarity  rather  than  speed  and  economy.  The 
program  Adamax3  was  written  for  demultiplexed  data.  The  operation 
could  be  speeded  up  by  operating  on  multiplexed  data. 

The  statement 

Call  Disc63 (NB,IRW, VAR ,NPRW) 

gives  a  read  (write)  from  (onto)  the  disc  if  IRW  is  0(1).  NB  is 
the  block  number  (32  words  per  block)  to  read  (write)  from  (onto) 
the  disc.  VAR  is  the  variable  location  in  core,  and  NPRW  is  the 
number  of  data  points  to  read  (write)  from  (onto)  the  disc. 

It  should  be  noted  that  subroutine  Upgrad  expects  the  filter 
vector  F(J,I)  to  be  stored  backwards  in  core  (i.e.,  F(L,I), 

F ( L- 1 , 1 ) _ F(1,I) . 


UU  U  OOUO  CJOO 


PKOGHAH  ADAMAX3 


part  of  the  program 


06  27  66 


this  is  ihf  main 


Ih*t>12-LF 

KK*iM/32 

Kl*(NpTS«i i/\H+l 

N2  *N“1 

Nj^atNpTS-!  )/324l 
JH«LF  ♦  } 

LF  lb  A  MU)  TIPLfc  OP  1J> 

DO  1U1  «sl.Kl 
DO  9  IaJ*N2 

9  CALL  DlSC6*< « l-lJ*Nl*<K-l>*KK*0*M < I-l>#5l2n>»5l2> 
Call  DISC!>j<N2*N1**K-1)*KK»Q'D*512> 
call  d i sc^  <  <n*i  1  *ni ♦  <k«i  i *kk  #o  * Y*5i2 ) 

Do  U  I *1. n2 
Do  13  J«l#5l2 
Lj«< 1-1 )*5i24J 
13  X(LJ*«D1  JUXlLJ* 

1FCK.EO«N1)104#1U5 
104  IM*NF'TS-<K1-1)  MM 


1q5  CONTINUE 

DO  1U0  rt*l»IM 
ISW2*H-1U0 

IFCK.GF.Ki-KCo)1&o»151 


*CO  IS  1  HE  CUTOFF  VA|  Ur  FOR  UPERATlNb  IN  MODE  1 
Last  POINT  FOR  MUuc  ,  nPERAT10N»tKl-HC0)*C5i2-LF> 

150  ISW2-1 

151  continue 

Call  UP(iMAn(X(H>#QPAn#F»N2» JM, Ti.SK* fc(H)»L#V2 #Ro* IS«2> 
MLF«M*LF 
Y ( HLF )  ■  Yj 

y(hlf*ii«y? 

E ( HLF )*UIMLF )«Yl 

100  6<MLF  *i >»D(MLF*1 »wY2 

call  op>c63<n*nimk-i )*kk*ll*i*e» jm>* jm> 

101  Call  DlSC63<  (N*l  I*N1  +  (K-1  )*KK4-LL*i  *  YC  JH)  .  JM) 


c 


o  o  o  o 


C 

C 


SUPHUlTHg  UPGPAUiVtn.F.N.LF.  y.!;k#e#l.Z,h0#  !SWg  J 
*HlS  SUbHOi'T  I NF  uPnATg«:  The  uradIennof  The  error 
A R p I  V A i.  uF  6ach  N6u  pATa  POINT, anD  CALCULATES  the 

old  «cf«thp  subroutine  then  outputs  The  Result  of 
FlLTfcR  WITH  The  iNPUT  DATa 


11 


13 

1« 

1* 

10 


DIMENSION  Y(^12*NI,QCL»N),F  (L,N)»E(LF  )  ,->  Rq  <N> 
lf«length  nf  fining  INTERVAL  ♦  i 
L*LfcNGTH  nF  FIL  >EP 
T«0  »u 
Zee  «U 

do  ii  jn,i 
Do  11  I«1#N 


MJ#l)«ti(j,iUX(i,|  >*E<  JUXILF-J  +  i ,  I  )*E(LF) 

Oo  10  I«1,W 


R0  *  1 J,Rfl  <  1 1  — x  ti • i !**?+y(LF» I >**2 
ski  *sk/Rij  <  i  >  € 

IFCISh?#LE.0»13#15 

DO  1«  J»1,L 

F  ( J,  I  )sF  O,  I  >aS^1*g(|  -J^i,  I ) 

°0  lo  .!■!#! 

T«Y*F  CJ* 1 )*X<J*LF«L* {  ) 

Z«Z+F ( J # 1 )*X( J+LF «L*1#  I ) 

return 

END 
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